Ensemble Theory for Force Networks in Hyperstatic Granular Matter 
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An ensemble approach for force networks in static granular packings is developed. The framework 
is based on the separation of packing and force scales, together with an a-priori fiat measure in the 
force phase space under the constraints that the contact forces are repulsive and balance on every 
particle. In this paper we will give a general formulation of this force network ensemble, and derive 
the general expression for the force distribution P{f). For small regular packings these probability 
densities are obtained in closed form, while for larger packings we present a systematic numerical 
analysis. Since technically the problem can be written as a non-invertible matrix problem (where the 
matrix is determined by the contact geometry), we study what happens if we perturb the packing 
matrix or replace it by a random matrix. The resulting P(/)'s differ significantly from those of 
normal packings, which touches upon the deep question of how network statistics is related to the 
underlying network structure. Overall, the ensemble formulation opens up a new perspective on 
force networks that is analytically accessible, and which may find applications beyond granular 
matter. 

PACS numbers: 45.70.-n, 46.65.+g, 83.80.Fg,05.40.-a 
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I. INTRODUCTION 

One of the most fascinating aspects of granular media 
is the organization of the interparticle contact forces into 
highly heterogeneous force networks flf. Direct evidence 
for these force networks mainly comes from numerical 
simulations 0, 0and experiments on packings of photo- 
elastic particles Ili_p. While the contact physics can be 
quite convoluted numerical studies have shown that 
qualitatively similar force networks occur in systems with 
much simplified contact laws |^ . It has nevertheless 
remained a great challenge to understand the emergence 
of these networks and their properties. 

Even though the spatial structure and anisotropics of 
the force network may be important [1, B S Hll j 
a more basic quantity, the probability density of contact 
forces P{f), has emerged as a key characterization of 
static granular matter Hi S Edi El) Il3 • Recently this 
quantity has also been studied for a wider range of ther- 
mal and athermal systems . Most of the attention 
so far has been focussed on the broad exponential-like 
tail of this distribution. Equally crucial is the generic 
change in qualitative behavior for small forces: P{f) ex- 
hibits a peak at some finite value of / for "jammed" 
systems which gives way to monotonic behavior above a 
glass transition [l^ ITtI . This hints at a possible connec- 
tion between jamming, glassy behavior and force network 
statistics, and underscores the paramount importance of 
developing a theoretical framework for the statistics and 
spatial organization of the forces |0| . 
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In this paper we study theoretical aspects of an en- 
semble approach that we recently introduced to describe 
these force networks This force network ensemble 

is based on the separation of packing and force scales 
that occurs in systems of hard particles: in most exper- 
iments, typical grain deformations range from 10^^ to 
10~^. The crucial observation is that these packings are 
usually hyperstatic, i.e., the amount of force components 
is substantially larger than the number of force balance 
constraints [l9j. This makes the problem "underdeter- 
mined" in the sense that there is no unique solution of 
the force network for a given packing configuration. For 
example. Fig. la shows two different force networks for a 
regular packing of 2D balls in a "snooker-triangle" . The 
ensemble is defined by assigning an equal a priori proba- 
bility to all force networks in which the net force on each 
particle is zero, for a given, fixed particle configuration. 
Since we want to describe non-cohesive particles, we then 
consider only those networks that have purely repulsive 
forces. As can be seen from Fig. 1, these simple rules 
indeed yield configurations that resemble realistic force 
networks, as well as a force distribution P{f) as typically 
observed in experiments and simulations. An important 
objective of this paper is to deepen our understanding of 
the force distribution, for this simplified but well-defined 
problem. 

Our ensemble approach is in the same spirit of the 
Edwards ensemble, in which an equal probability for 
all "blocked" or "jammed" configurations is postulated 
pol l2l| . This Edwards ensemble does not only average 
over forces, but also over all possible packing configura- 
tions, which makes the problem difficult to track theo- 
retically. We therefore propose to exploit the separation 
of length scales that occurs for hard particles, by fixing 
the packing geometry (macroscopic scale) and allowing 
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FIG. 1: (a-b) Two different mechanically stable force con- 
figurations for a "snooker-triangle" packing of 210 balls; the 
thickness of the lines is proportional to the contact force. The 
"force network ensemble" samples all possible force configu- 
rations for a given contact network with an equal probability, 
(c) After sampling many force configurations, this yields the 
following distribution of interparticle forces P{,f)- 



for force fluctuations (microscopic scale). Besides prac- 
tical advantages, the conceptual gain of separating the 
contact geometry from the forces is that we can start 
to disentangle the separate roles of contact and stress 
anisotropies |3, HO) uM ■ Interestingly, the idea to restrict 
the Edwards ensemble to fixed packing geometry has also 
been proposed recently by Bouchaud in the context of 
extremely weak tappin g [ 221 . and was also employed in 
recent simulations |23l |2J| . Note that this force ensem- 
ble incorporates the local force balance equations on all 
particles and therefore it is fundamentally different from 
recent entropy-based models for force statistics |2^ . 
In these studies one postulates an entropy functional in 
terms of the single force distribution P(/), without in- 
cluding the intricately coupled force balance equations 
and resulting force correlations. 

From a more general point of view, the ensemble pro- 
vides a challenging statistical physical problem of rather 
broad interest, that of sampling the solution space of a 
set of underdetermined equations and constraints. For 
example, the problem is mathematically very similar to 
the so-called flux balance analysis that is used to unravel 
metabolic networks in biological systems l28l|. Here 
the reaction fluxes are underdetermined and play a role 
analogous to the forces discussed here. In contrast to the 
forces, however, these fluxes typically display power-law 
distributions |2a |. This touches upon the deep question 
of what kind of statistics emerges when balancing scalars 
on a network of a given structure 0, , and shows that 



the nature of the set of balance equations has a strong 
influence on the resulting statistics. 

The aim of this paper is to explore the "phase space of 
force networks" and to unravel how this gives rise to the 
robust characteristics of the force distribution P(/). We 
will initially focus on regular packings which are highly 
coordinated and therefore far from the isostatic limit. 
The advantage of these packings is, however, that the 
underlying physics is more transparent and that small 
regular packings can be resolved analytically. In addi- 
tion, their force distributions are quite comparable to 
those found in numerical explorations of the ensemble 
for amorphous packings presented elsewhere 0, |33| • We 
will also study the ensemble on generalized networks, for 
which the force distributions rapidly loose their similarity 
to those of real packings. 

After defining the ensemble in more detail, the paper 
consists of four parts. In Sec.m we study the force en- 
semble for spherical, frictionless particles in regular tri- 
angular "snooker" packings. We discuss how these force 
distributions are related to geometric aspects of the high- 
dimensional phase space. In Sec. Illll we provide a formal 
mathematical description of the ensemble and derive the 
explicit form of P{f), Eq. This expression contains 
coefficients that depend on the packing geometry, and 
which we have been able to compute for several small sys- 
tems. These exact P{f) already exihibit the features that 
are relevant for larger, more realistic packings, and will 
be presented in Sec. II VI Due to the linearity of the equa- 
tions of force balance, the problem can be further gen- 
eralized by considering perturbations of the packing ma- 
trix and random matrices, which are presented in Sec.lVl 
This probes which ingredients are essential for obtaining 
realistic P(/)'s. The paper closes with a discussion of 
the strengths and weaknesses of our approach and indi- 
cates some open issues and other problems that can be 
addressed with the ensemble. 



A. Definition of the force network ensemble 

We will now introduce the main aspects of the ensem- 
ble approach. Even though our approach is perfectly 
suited to include frictional forces [ill l2^ l24l | . for simplic- 
ity we will restrict ourselves to packings of N frictionless 
spheres of radii Ri with centers . We denote the inter- 
particle force on particle i due to its contact with particle 
j by iij. There are zN/2 contact forces in such packings 
{z being the average contact number), and for purely 
repulsive central forces we can write — fijVij /\rij\, 
where all fij {— fji) are positive scalars. For a fixed 
contact topology in d dimensions, we are thus left with 
dN unknown positions and zN/2 unknown forces fij. 
Note that the number of unknown forces is not precisely, 
but close to, zN/2 if boundary forces are present. 

These degrees of freedom satisfy the conditions of me- 
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chanical equilibrium, 
dN eqs.: ^ ftj = 0, where = - rj, (1) 

and once a force law F is given, the forces are explicit 
functions of the particle locations: 

zN/2 eqs.: /„■ = F{r,j ■,R^,Rj) . (2) 

The contact number z is a crucial quantity. As has 
been argued before 0, |^ , even though packings of 
infinitely hard frictionless particles have z = 2d and are 
thus isostatic, for particles of finite hardness, packings are 
typically hyperstatic with z > 2d. In this paper we focus 
on hyperstatic packings, but before doing so, we wish 
to point out an important subtlety. In recent numerical 
work, it was shown that z approaches the isostatic limit 
for vanishing pressures (hence vanishing deformations) of 
the particles, and that the (un)jamming transition here 
is similar to a phase transition, with power-law scaling 
of the relevant quantities and the occurrence of a large, 
possibly diverging lengthscale [s^. Therefore, the precise 
value of z may be important, since it reflects the distance 
to the jamming phase transition; this may bear on the 
interpretation of our results. It is worth pointing out that 
for frictional packings, even in the limit of infinitely hard 
particles, z stays away from the isostatic limit .34]. 
Hyperstatic packings are therefore important, and our 
work, even though it focusses on frictionless packings, 
may also be seen in this light. 

Returning to the force network ensemble, in the regime 
where particles are hard but not infinitely hard, varia- 
tions of the force of order (/) result in minute variations 
of rij. Hence, Eqs. (Q) and Q can effectively be con- 
sidered separated, and the essential physics is then given 
by the force balance constraints Eqs. with fixed r^. 
In this interpretation there are more degrees of freedom 
{zN/2) than constraints (2A^), leading to an ensemble of 
force networks for a fixed contact geometry. 

This ensemble for a fixed contact geometry is then con- 
structed as follows, (i) Assume an a-priori flat measure 
in the force phase space {/}. (ii) Impose the 2N linear 
constraints given by the mechanical equilibrium Eqs. 
(in) Consider repulsive forces only, i.e., V/y > 0. (iv) 
Set an overall force scale by applying a fixed pressure or 
fixed boundary forces, similar to energy or particle num- 
ber constraints in the usual thermodynamic ensembles. 

We are thus considering the phase space defined by 
the force balance Eqs. (Q, the condition that all /'s are 
positive, and a "pressure" constraint J2kfk = Ptot- For 
notational convience, we indicate the forces by a single 
index k throughout the remainder of paper. Since all 
equations are linear, the problem can be formulated as 

Af^b and V /fc>0, (3) 

where the fixed matrix A is determined by the pack- 
ing geometry, / = {fi, f2, ■ ■ ■ JzN/2), and b = 
(0,0,0,---,0,Ftot). 
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FIG. 2: Three monodisperse frictionless spheres in a snooker 
triangle. This system has 9 unknown forces: 6 boundary 
forces (/i to /e) and 3 interparticle forces (/t, fs, /g). 



II. REGULAR PACKINGS: BALLS IN A 
SNOOKER TRIANGLE 

In the introduction we have seen that our ensemble 
approach for a snooker packing of 210 particles repro- 
duces a force distribution that is very similar to those 
obtained in experiments and simulations. To understand 
how this shape of P{f) comes about, we now work out the 
force network ensemble for small systems of crystalline 
(monodisperse) packings. We first study the packing of 
3 balls shown in Fig. [51 for which we explicitly construct 
the phase space of force networks. As this system is very 
small, the force distribution deviates considerably from 
distributions observed in large systems. It nevertheless 
provides a very instructive example. We then present a 
numerical analysis of how P{f) evolves as a function of 
system size for snooker packings. Remarkably, a packing 
of 6 balls is already sufficiently large to obtain the char- 
acteristic peak in P{f)- We therefore address general 
physical aspects by elaborating on this system. 



A. Three balls 

In the system of 3 balls depicted in Fig. [3 we encounter 
9 unknown forces: 6 boundary forces and 3 interparticle 
forces. These forces have to balance on each particle in 
both the X and y directions, which constitute 2x3 = 

6 linear constraints. In addition, we impose an overall 
pressure by keeping the total force on a boundary at a 
fixed value: for example we fix fi + f2 = 2. Interestingly, 
one can show that such a boundary or pressure constraint 
is equivalent to keeping the sum over all forces at a fixed 
value: in appendix^we demonstrate that keeping fj 
at a constant value is equivalent to a constant pressure, 
also for irregular packings. 

Together with the pressure constraint, there are thus 

7 linear equations to determine the 9 unknown forces 
/ = (/i; , ■ • • , /g), and hence, there is a two-dimensional 
space of solutions. This space does not contain the origin 
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FIG. 3: (a) The 2D phase space of the 3-baU problem can 
be defined using three simple independent solutions of the 
problem, (b) The first solution /a has /i = /4 = 2, f5=fe = 
1 and /g = \/3 and /2,3,7,8 = 0; the other solutions fs and 
fc follow from the threefold symmetry of the packing. 

of the force space, for which all fj = 0, due to the in- 
homogeneous pressure constraint. As a consequence one 
requires three vectors to characterize the two-dimensional 
space: two basis vectors and a vector defining the loca- 
tion of the plane with respect to the origin. Using linear 
algebra one can construct these three vectors from three 
linearly independent force network solutions /a, /s, /d 
which allows to express the general solution as 




FIG. 4: Two dimensional cut through the phase-space 
spanned by the nine forces of the 3-ball problem, (a) The 
borders of the triangle are the lines where one of the inter- 
particle forces changes sign; the shaded area represents the 
probability to find a configuration between and /t + Sfj. 
(b) The corresponding force distribution P{f7). 




FIG. 5: Two dimensional cut through the phase-space 
spanned by the nine forces of the 3-ball problem, showing 
how boundary forces are distributed, (a) The borders of 
the hexagon are the lines where one of the boundary forces 
changes sign; the shaded areas represents the probability for 
a certain /i. (b) The corresponding force distribution -P(/i). 



/ = caIa + cbJb + {1-ca- CB)fc ■ (4) 

An intuitive picture of this equation is provided in 
Fig. the 2-dimensional plane can be defined from 
three solutions (very much like a line can be defined by 
two points). However, the constraints that all fj > 0, 
provide serious limitations on the allowed values of ca 
and cb- As will be shown below, only a small convex sub- 
set of the the two-dimensional solution space represents 
force networks consisting of strickly repulsive forces. 

Using the solutions of Fig. 13 to construct the phase 
space, we obtain the triangle depicted in Fig. In this 
picture, the three solutions are the corners of the triangle. 
For example, the right corner represents the first solution 
in Fig. 01 /a for which {f7,hJj>) = (0,0, \/3), whereas 
the left corner corresponds to /s that has {fjjfsifd) — 
(0, a/S, 0). A superposition of these two vectors is still a 
solution of our linear problem, and since in both cases 
/7 = 0, the base of the triangle is a line where /y = 0. 
The upper corner represents fc, for which f^ attains 
its maximum value of -n/S- Therefore, the dashed line 
is a projection of the fr axis onto this 2D space of so- 
lutions. This implies that the space below the triangle 
corresponds to a region where /y < 0, which is forbidden 
for repulsive particles. Applying the same argument for 



fs and /g, one realizes that only the area inside the tri- 
angle is allowed. As we mentioned in the introduction, 
the ensemble assumes an equal a priori force probabil- 
ity, which makes each point in the triangle equally likely 
(due to the linearity of the force balance restrictions). 
Therefore, the probability to have a solution between /y 
and fr + Sfj is simply represented by the shaded area 
in Fig. This "volume" decreases linearly as /v ap- 
proaches its maximum value, so that the distribution of 
fj simply becomes P{fr) = |(^/3 - /7)e(\/3 - fj) - see 
Fig.Et). 

The combination x&{x), where &{x) is the Heaviside 
step function, will occur in most P{f) thoughout this 
paper. Therefore we introduce 

T{x) EE xe{x) . (5) 

The distribution of the boundary forces (/i to /g) can 
be found in a similar manner. Checking the three in- 
dependent solutions, one finds that /i = at the left 
corner, /i = 1 at the upper corner, and /i = 2 at the 
right corner of the "phase space triangle" . From the geo- 
metric construction in Fig.lSJt, it is easy to find the /i = 
line, and the projection of the fi axis is indicated by the 
arrow. Due to symmetry there are of course six such 
borders (/i = to /g = 0), and all boundary forces are 
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FIG. 6: P{f) for bulk forces in "snooker" packings of increas- 
ing sizes. 



positive inside the hexagon. So, the solutions for which 
all forces are positive lie within the triangle. Considering 
the shaded area in Fig. we obtain the distribution of 
boundary forces P(/i) = T{2 - fi) - 2T{1 - fi), which 
is shown in Fig. (SJ). We thus find that there is a qualita- 
tive difference between the boundary forces (/i, • • • , /e) 
and the interparticle forces (/y, fs, /g). Interestingly this 
is also the case for larger systems and is consistent with 
earlier work on statistics of wall forces 36, 37] . 

Although this 3-ball system provides a very nice illus- 
tration of how to obtain P{f) from all possible force con- 
figurations, it is not complex enough to reproduce non- 
monotonic P{f)- In fact, the problem discussed above 
is equivalent to partitioning a conserved energy into 3 
positive parts. In our case, the conserved quantity is the 
total force and the 3 parts are the coefficients ca, cb and 
(1 — ca — cb)- In the thermodynamic limit, the problem of 
partitioning, e.g., an energy E^ot simply yields the Boltz- 
mann distribution of energies Ei of subsystems; also for 
finite systems these distributions are always monotoni- 
cally decreasing - see appendix In Sec. Ill ('I we show 
that the problem of 6 balls already has enough complex- 
ity that it leads to non-monotonic behavior of P{f)- 



B. Numerical analysis of larger systems 

To compute P{f) for larger packings, we have applied 
a simulated annealin g pr ocedure |35| . As was shown 
in our previous work [TJ this scheme can also be used 
for irregular packings. Starting from an ensemble of 
random initial force configurations taken from an arbi- 
trary distribution with (/) = 1 and fj > 0, we se- 
lect a random bond j and add a random force A/, so 
that fjiji) — fj{o) + A/, in which the symbols n and 
o denote the new and old force respectively. The ran- 
dom change from o to n is accepted with a probability 
given by the conventional Metropolis rule p{o n) = 
min {l,e (/j (n)) exp [- (H(n) - H(o)) /T]), in which Ti. is 
a penalty function whose degenerate ground states are 
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FIG. 7: Boundary forces for "snooker" packings of increasing 
sizes. 



solutions of Eq. 

n{f}^[Af-by. (6) 

For large packings {N > 500) it is computationally much 
more efficient to always satisfy (/) = 1 by selecting two 
bonds (j 7^ k) at random and using fj(n) = fj{o) + Af 
and fk{n) — fk{o) — Af as the update scheme, so that 
the pressure constraint can be left out of the penalty 
function. By slowly taking the limit of T ^ we sam- 
ple all mechanically stable force configurations for which 

— > 0. We have carefully checked that results do not 
depend on the initial configurations and details of the 
annealing scheme. In section Hvl we will show that this 
scheme perfectly reproduces analytic results for small 
regular packings. 

The two force networks shown in Fig.^are typical solu- 
tions / obtained by this procedure. The resulting distri- 
butions of interparticle forces are presented in Fig. El for 
packings of increasing number of balls; boundary forces 
will be discussed seperately and are not included in these 
P(/). Note that all P(/)'s display a peak for small /, 
which is typical for jammed systems llg. The fact that 
the probability for vanishing interparicle forces remains 
finite is in agreement with most numerical and experi- 
mental observations; only a few studies report power-law 
behavior for small forces 0,^3- ^'^^ large packings, this 
peak rapidly converges to its asymptotic limit. The tail 
of P{f) broadens with system size, and will be discussed 
in more detail in Sec. I VII 

In Fig.[7|we show the probability distributions P(/„aii) 
for the forces /wall between sidewalls and balls for the 
regular packings of increasing size. As has been discussed 
at length in [3a. l37l | , these distributions differ from the 
probability distributions of bulk forces. In this particular 
case this is easy to see: the boundary force /wall has to 
balance the force of the two balls in the next layer with 
which it makes contact (excluding the corner balls). Even 
though each of these forces has a finite probability to be 
vanishingly small, the probability that both these forces 
are small has not, hence P(/waii) for /waii 0. 
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C. P{,f) and phase space geometry 

Here we will discuss some geometrical aspects of the 
set of allowed force configmations. Consider the zN/2 
dimensional force phase space spanned by the fj. Since 
all 2N force balance equations iQJ are linear in the forces, 
the allowed solutions lie on a hyperplane of dimension 
{zN/2 — 27V). (Note that the overall pressure constraint 
introduces an additional constraint, lowering the dimen- 
sion by 1.) Furthermore, since we consider repulsive 
forces only, this plane is restricted to the "positive" hy- 
perquadrant where all fj > (see Fig. IHJ. Therefore, 
the allowed force-configurations form a (hyper)-polygon 
whose facets are given by the conditions that some force 
fj becomes 0. Under our assumption of a "flat measure", 
all points on this polygon correspond to valid force net- 
works with equal probability. 

A number of basic properties of this solution space can 
now readily be deduced. Trivially, the solution space is 
convex: due to the linearity of the equations, the points 
on a straight line connecting two admissible solutions are 
admissible solutions themselves, as was also pointed out 
in I23. l33 |. Although this is not immediately ob- 
vious in low dimensions, for higher dimensional bodies 
the overwhelming part of the "measure" is concentrated 
near the boundary (think of a high dimensional sphere, 
where almost all volume is in the "shell" close to sur- 
face). Near the boundaries, one or more forces tend to 
zero, and this is consistent with the fact that in typi- 
cal force networks a finite fraction of the forces are close 
to zero (since P{f | 0) 7^ 0). More homogeneous force 
networks, for which all forces are around some average 
value, correspond to points in the phase space that are 
sufficiently far away from the boundary. While such con- 
figurations are perfectly allowed within our framework, 
and are easy to construct by considering a suitable lin- 
ear combination of "ordinary" force- networks, they only 
occur with vanishingly small probability in the limit of 
large N, and are thus extremely unlikely to be seen in 
"unguided numerics" or experiments. 

Even though we have not worked this out in detail, we 
expect that some more general properties of the force net- 
works could be related to geometrical properties of (ran- 
dom) hyperpolygons. As one simple example consider 
the following. For two forces, say fi and fj, to become 
zero simultaneously, the facets i and j have to touch; in 
general this may not be possible geometrically, so that an 
intruiging issue concerning correlations between distant 
forces arises. 

Another issue that may have a relatively simple inter- 
pretation in the polygon language is the peaked appear- 
ance of P{f )- We suggest the following intuitive picture, 
based on a consideration why the slope dP{f)/df can 
be expected to be positive for small forces. For very 
small systems, like the case of 3 balls discussed in section 
III Al this is not true. This immediately follows from the 
shape of the allowed phase space polygon. As shown in 
Fig. ^ this is a triangle where the angles between the 




FIG. 8: Schematic representation of the phase space of al- 
lowed force configurations. Each fij defines a direction in the 
zN/2 dimensional force space. By imposing the linear condi- 
tions of mechanical equilibrium, this space is restricted to a 
"hyperplane" of lower dimensionality. The physically allowed 
region is a (hyper)polygon is bounded by the requirement that 
all fij > 0. 



bounding edges were acute. When we move away from 
a / = boundary, the phase space volume decreases so 
that dP/df < 0. If we go to larger systems, however, the 
number of facets bounding the space {— zN/2) becomes 
much larger than the dimension D (= zN/2 — dN — 1) of 
the polygon. Hence, we expect that the "angles" between 
bounding facets will typically become obtuse, which will 
make the phase space increase when increasing /. This 
indicates that dP/ df is typicallypositive for small forces, 
so that P{f) displays a peak |39| . 

Six halls - Let us provide another perspective on the 
phase space geometry by discussing the problem of 6 
balls, which is the smallest snooker packing displaying a 
non-monotonic P{f)- For the 6 balls there are 18 forces, 
which are constrained by 2 x 6 -I- 1 = 13 equations, so the 
space of solutions is a 5D hyperplane. If we try to con- 
struct the phase space like we did for the 3 balls, we now 
require 6 independent solutions / that obey force balance 
on each particle. Again, there exist simple solutions of 
linearly propagating force lines - see Fig. |^. However, 
there are only 3 such solutions, so we also require nontriv- 
ial solutions where forces "scatter" at a certain particle. 
For example, we can take 3 solutions of the type shown 
in Fig. EId. 

The presence of these nontrivial solutions changes the 
phase space in a a fundamental manner. A given force 
can now take a certain value in many different ways, by 
different linear combinations of elementary "modes" . In 
other words, a force can no longer be associated to a sin- 
gle mode of the force network, like it was the case for the 
three interparticle forces in Fig. |2| As a consequence, 
the problem has become much more intricate than sim- 
ply partitioning the total force into positive amplitudes 
(which, for large systems, would lead to a simple expo- 
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(a) 




where the nonzero aij are projection factors between — 1 
and +1. There are dN such equations, which we label 
as z = 2, 3, • • • , dN + 1. To keep the overall pressure at 
a fixed value we impose fj — F, which for notational 
convience we write as 





FIG. 9: Two different types of solutions of force equilibrium 
for 6 balls in a snooker-triangle. 



nential distribution, see appendix EJ- Instead one finds 
nontrivial force distributions, for which we derive analyt- 
ical expressions in the following section. Indeed, for all 
investigated packings, we observe non-monotonic P{f) 
whenever scatter-solutions occur. 



III. GENERAL FORMULATION FOR 
ARBITRARY PACKING GEOMETRY 

In this section we show how statistical averages can be 
computed analytically within the force network ensemble, 
for arbitrary packings. We present a systematic way to 
evaluate the complicated high-dimensional integrals as a 
sum over contributions of the following structure: 



zN/2 

aijfj — F , with all aij = 1 . (9) 

We thus encounter a matrix problem Af = b, where the 
the components of A. 
Imposing the various constraints and assuming an 
equal a priori probability in the force space defined by 
/ = (/i, /2, • ■ • , fzN/2), we obtain the joint probabihty 
density 
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which is normalized by the phase space volume 
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P{f) = Ec;,rMl-&A./)^"'-'^^e(l-6A/) 



= 5]c;,rM^(l-&A/)] 



D-l-qx 



Since we consider repulsive forces, / df represents an 
integral over all forces in the hyperquadrant where all 
fj > 0. With this measure, we can now compute the 
(7) single force distribution P{fj) as 



where D is the dimension of the phase space, and the co- 
efficients bx, c\ and qx depend in a nontrivial way on the 
particle packing; for most A, we find that qx = Q. The 
function T was defined in Eq. jSj) ; note that the contribu- 
tions [T(l — bfyf~^ ^ g~{D~i)bf jj-^ ^j^g thermodynamic 
limit. For the reader who is interested in the results but 
not in the details of the derivation, we summarize exact 
P{f) for small regular packings in Sec. IIVI 



A. Mathematical definition of the ensemble 

The phase space of force networks is defined by the lin- 
ear constraints of force balance, an inhomogeneous lin- 
ear constraint to fix the pressure, and the requirement 
that all forces are non-negative. If we now indicate each 
contact force by an index j, we can express mechanical 
equilibrium as 



zN/2 



(8) 



n 



dfk 
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(12) 



which in principle can be different for each fj] for exam- 
ple see the boundary forces within the snooker-triangles 
(Sec.^. In practice, it turns out that P(/j) for different 
interparticle forces shows only little variation. 

The fact that we only integrate over the hyperquadrant 
where all fj > makes it difficult to evaluate the integrals 
explicitly: each integration of the 6 function gives rise to 
a Heaviside 9 function to keep track of the boundaries 
of the phase space. To avoid this problem we represent 
the 6 functions as Fourier integrals 
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dsi_ 
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(13) 



which has the advantage that the fj only occur in an 
exponential way and they are easily integrated out. If 
we now write s = ^^(si, ■ ■ • , Sm), where m = dN + 1, the 
partition function Q becomes 
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where the factor e^''^^ arises due to the inhomogeneous 
pressure constraint We furthermore added cut-off 

factors e~'^J so that the integrations over the fj are defi- 
nite; at the final stage we take the limit tj 0. The rows 
of the matrix A correspond to the constraint variables Si 
and the columns correspond to the denominators origi- 
nating from the fj integrals. From now on we indicate 
the dimensions of the matrix hy m — dN + 1 (number of 
rows) and n = zN/2 (number of columns). 

All integration variables Si run from — oo to oo, so we 
can evaluate them as contour integrations in the com- 
plex plane. The integrand is a product of denominators, 
and each Si occurs in as many denominators as there are 
forces acting on a certain particle. In the absence of grav- 
ity, each mechanically stable particle should at least have 
3 contacts. This makes the integration over the Si con- 
verging at infinity and allows to close the contour either 
through the upper half plane or through the lower half 
plane. An exception is the si integration, which has to 
be closed through the upper plane since F > 0. 

Let us first integrate out Sm- Each denominator that 
has Umj 7^ gives rise to a pole at 



Smij) 



E 



(15) 



The residue is obtained by substiting this pole in the 
remaining n — 1 denominators of Eq. H14(l . Note the im- 
portance of the ej to make the integration definite. It 
is easily seen that this substitution leads to a renormal- 
ized matrix A* of m — 1 rows (constraint variables) and 
n—l columns (denominators), and to renormalized e^/^j 
as well. However, the key observation is that the remain- 
ing integrals are of the same type as Eq. H14|) . We thus 
find a recursion relation 



^mn(-4) ^ ^ ^ ^ni—l.n—l{Aj) 



(16) 



where the sum extends over all encircled poles. The sym- 
bol ± indicates that the contribution is positive or neg- 
ative depending on whether the integral has been closed 
through the upper (+) or lower (— ) half plane. The 
renormalization to A* is different for each pole, so each 
term has to be followed independently. At each integra- 
tion the number of contributions therefore grows rapidly, 
since each new pole gives rise to a new "branch" of the 
recursion Eq. (|16|l . The exponential increase of the num- 
ber of branches with the size of the tree forms a severe 



limitation on the solutions for larger systems. At the 
final stage, we have to compute = ^i.d+i by 

integrating over si: 
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where D is the dimensionality of the phase space. The aij 
and ej appearing in this equation are obtained from suc- 
cessive renormalization each time a pole is substituted. 

So, the calculation of Q involves a tree-like structure 
where the branching rate is equal to the number of en- 
circled poles. Using relation Eqs. (|16|l and H17|l one can 
compute the contribution of each individual branch, us- 
ing a recursive scheme. The fact that fi scales as F to the 
power D is not surprising: F is the only force scale for 
the D dimensional phase space, and in fact, the behavior 
F^ is obtained immediately from a trivial rescaling of 
Eq. I|ll|l . However, in the following paragraphs we show 
how the analysis presented above can be extended to the 
nontrivial calculation of the force distribution P{f). 



B. Calculation of P{f) 

Comparing Eqs. 1)11(1 and (|12|) . we notice that the ex- 
pression for P{fj) is the same as that for Vl without the 
integration over fj ; without loss of generality we will con- 
sider P{fi). As a result, the expression for P{fi) contains 
one less dominator than Eq. (|14|l and instead there will 
be an additional exponential factor, i.e. 



(18) 



Following the same integration strategy as for fi, we again 
obtain a recursion of the type 

p„,„(/i) = ±;^— p™_i,„_i(/i) , (19) 



where for clarity in notation we left out the explicit de- 
pendence on the (renormalized) matrix A. After succes- 
sive substitution of the poles, the final integration over 
si becomes 
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Each branch of the tree gives a contribution of this type 
and together they accumulate to the result of Eq. O 
with qx = 0. The coefficients 5a are thus simply the 
ail/ F that remain after successive renormalization of the 
matrix A. We will demonstrate that, fortunately, the 
final result contains only a few different b\, at least for 
small packing geometries. 

In the final integration of Eq. H20|l. we implicitly as- 
sumed that all aij appearing in the denominators are 
not equal to zero. They may become negative, provided 
that the associated small ej is also negative so that the 
pole is still in the upper half plane and the integration 
remains finite. Naively one would expect that it very 
unlikely that some aij = 0, since it corresponds to an 
accidental coincidence of two poles. However, for regu- 
lar structures like the snooker-packings it is a frequently 
occuring phenomenon. The double poles are responsible 
for the cases q\ ^ 0. We have adapted the algorithm 
such that it can deal with an arbitrary multiplicity of 
the poles. In some cases, these multiple poles alter the 
general result for P{f) with additional contributions of 
the type 
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These contributions can be recognized as the q\th deriva- 
tives of the general result, corresponding to the coinci- 
dence oi q\ + l poles. We expect, however, that multiple 
poles will never occur for disordered packings. 



IV. EXACT RESULTS FOR SMALL 
CRYSTALLINE PACKINGS 

We now present a number of exact P{f) for small crys- 
talline packing geometries. In particular, we have worked 
out the problem of 6 balls in a snooker-triangle, triangu- 
lar 2D packings with periodic boundary conditions, as 
well as a small 3D FCC packing with periodic bound- 
ary conditions. Following the algorithm described in the 
previous section, we have been able to obtain the coeffi- 
cients b\ and c\ appearing in Eq. {Tj) for these systems. 
For notational convenience, we express the results in the 
dimensionless force x — f /F. All is in perfect agreement 
with our numerical simulations. 

The intricate combinatorics has been performed using 
a computer program. As mentioned the number of con- 
tributions grows exponentially with the size of the tree, 
since the branching rate is of order of 2 per elimination 
step. Even worse is the fact that the different signs of 
the contributions lead to large cancellations. The results 
given below for small systems are the result of many more 
terms in the tree. This makes the algorithm numerically 
unstable for larger systems. 
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FIG. 10: (a) All solutions of the 2x2 periodic arrangement 
can be described as a superposition of linearly propagating 
force- lines, (b) The corresponding monotonic P{f)- 



A. 2D triangular packings with periodic boundaries 

Four balls - The smallest interesting 2D triangular 
packing with periodic boundary conditions is the 2x2 
packing of 4 balls. It has 3 x 4 = 12 unknown forces and 
2x4 = 8 equations expressing mechanical equilibrium. 
Due to the periodic boundaries, however, two of these 
equations are actually dependent. Hence there are only 
6 independent equations and together with the overall 
pressure constraint this results into a, D = 12— (6+1) = 5 
dimensional phase space. 

In terms of the dimensionless variable x — f/F, we 
obtained the following result for this system: 



P{x) = 10T^(1 - 2x) 



(22) 



Taking F = 12 so that (/) — 1, we plotted this dis- 
tribution in Fig. llUb . It is a monotonically decreasing 
function that allows a maximum force of Xmax = 5, i.e., 
/max = ■^F. This maximum force is achieved for a sim- 
ple "propagating" solution shown in Fig. IlUb : the total 
force F is shared between two nonzero forces only (note 
the similarity to the solutions shown in Fig. |3| for the 
packing of three balls). Due to the symmetry of the 
problem there are six such trivial solutions, which are 
in fact sufficient to define the whole 5D phase space of 
force networks. The 2x2 problem is therefore equivalent 
to partitioning the total force into six nonnegative "am- 
plitudes", just as was the case for the three balls in the 
snooker-triangle. Indeed, Eq. H22|l is of the same form as 
Eq. ljA6|l in appendix El 

Nine balls - For the 3x3 packing of 9 balls there are 
3x9 = 27 unknown forces that are constrained by 2 x 9 — 
2=16 independent equations of mechanical equilibrium. 
Fixing the overall pressure, one is left with a D = 27 — 
(16-1-1) = 10 dimensional phase space. This space can not 
be reconstructed from the trivial propagating solutions, 
of which there are only 9. Again, the presence of the 
"scatter" solutions such as the one shown in Fig. Illh . 
results into a non-monotonic P{x): 



P{x) = 40 



r9(l-3x)-^T9(l-9x) 



(23) 
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FIG. 11: (a) The system of 3 x 3 balls allows for nontrivial 
"scatter" solutions, (b) The corresponding is therefore 

non-monotonic. The solid curve is Eq. 1281 the crosses are 
obtained from numerics as described in section ITl Bl 



Taking = 27 so that (/) = 1, we plotted P{f) as a solid 
curve in Fig. Illb : the crosses indicate the distribution 
obtained by the same numerical method that was used for 
the snooker-triangles in Sec. The perfect agreement 
illustrates the accuracy of our numerical method. 



B. 3D FCC packing with periodic boundaries 

To illustrate that our ensemble can be applied to three 
dimensional packings just as well, we have computed 
P{f) in the conventional FCC unit cell, with periodic 
boundary conditions. This is a system of 4 balls, since 
the FCC unit cell contains 8 particles at corners (each 
counting for 1/8) and 6 particles at the faces (counting 
for 1/2). The coordination number of the FCC packing 
is z = 12, so there are zN /2 = 24 forces in this system. 
We now have to respect force balance in three dimen- 
sions, i.e. 3x4 = 12 equations, of which, due to periodic 
boundary conditions, only 9 turn out to be independent. 
Together with the pressure constraint, there are thus 10 
equations to constrain 24 forces, and hence the problem 
has a 14 dimensional space of solutions. 

The resulting P(x) turns out to be 



T"(l-2x)- Ari3(i_6x) 



— Ti3(l - 8a;) - 21xT^'^{l - 6x) 



,(24) 



Fig. 1 121 shows that this force distribution has the same 
typical features as those obtained for two-dimensional 
packings. It's a non-monotonic function, which can again 
be related to the existence of "scatter" solutions. There 
are 15 independent solutions to fix the 14D phase space 
of force networks, 12 of which are linearly propagating 
"trivial" solutions (two for each lattice direction). The 
other three are again "scatter" solutions. One of these is 
shown in Fig. E| 
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FIG. 12: (a) One of the "scatter" solutions for the FCC unit 
cell with periodic boundary conditions. The black spheres be- 
long to this unit cell; the grey spheres belong to neighbouring 
cells. All forces have the same magnitude; those within this 
unit cell are drawn as thick solid lines; the others are drawn 
as thick dashed lines, (b) The corresponding non-monotonic 
P{f), from Eq. 1241 (solid curve) and from numerics as de- 
scribed in section Hi Bl (crosses) . 



We already showed that one has to distinguish between 
the interparticle forces and the particle-wall forces, which 
obey qualitatively different statistics. Upon closer in- 
spection, however, one notices that there are also two 
different types of interparticle force: the 6 closest to the 
boundary (type /) and the 3 closest to the center (type 
//). We find that 
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where a = 2 {l + Vs). 

The numerical results shown in Fig. were obtained 
without discriminating between type / and type //. This 
is allowed since even though Pi{x) and Pii{x) are not 
identical, their shapes are very similar. Comparing the 
data with ^Pj{x) + ^Pii{x) gives again an excellent 
agreement between the theoretical result and the numer- 
ical result shown in Fig. El The factors 2/3 and 1/3 
appear because there are 6 forces of type / and 3 forces 
of type //. 

Finally, let us discuss the statistics for the boundary 
forces as shown in Fig. [3 Also in this case there are 
two different types of boundary forces, namely 6 at the 
corners (c) and 3 in the middle (m) of each boundary. 
We find that 



C. 6 balls in a snooker-triangle 

We now provide the exact force distributions for the 6 
balls in a snooker-triangle, which we discussed in Sec.lTTI 
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54(7 + 4^3) 

[-T^{l-bx) + 8bxT^{l-bx) + T-^{l-2bx)] , (28) 

where b — 3 + \/3. The hnear combination |Pc(a;) + 
^Pm{x) fits the boundary force distributions as shown in 
Fig. H extremely weU (not shown). 

V. BEYOND PACKINGS 

In the preceding sections we have extensively studied 
the force distributions emerging in the ensemble of force 
networks, for a variety of crystalline packings. The var- 
ious P{f) are non-monotonic and display only marginal 
differences. As we demonstrated in Ref. (ll|, the same 
qualitative behavior is observed for irregular packings. 
Even though the packing matrices differ substantially in 
these cases, the resulting P{f) is extremely robust. This 
raises the question of which are the essential ingredients 
to obtain a typical force distribution. In other words, 
what properties of the packing matrix A determine the 
shape of P(/)? 

All packing matrices consist of a large number of zeros, 
except for a few elements per row that are projection fac- 
tors between — 1 and 1 . Such a matrix has some features 
of a random matrix, but it implicitly contains the en- 
tire spatial structure of the system. To see whether this 
spatial structure is crucial for the typical shape of P{f), 
we now study true random matrices, which no longer 
represent a physical packing of particles. Of course, we 
still extend the matrix by the normalization constraint 
fj = F and demand that all fj>0. 

We find that such random matrices yield P{f) whose 
decay is described by a product of Gaussian and expo- 
nential tails. However, all these distributions are mono- 
tonically decreasing and thus lack the typical peak, even 
when considering "sparse" random matrices. We then 
try the opposite approach, where we start from a phys- 
ical packing matrix and then slowly introduce random- 
ness. In contrast to the striking robustness of P{f) for 
real packings, the force distribution is very sensitive even 
to small perturbations away from the physical matrix. 



A. Random matrices 

1. Infinite Gaussian random matrices 

We start out the random matrix approach by generat- 
ing all elements in Eq. JS)) from a Gaussian distribu- 
tion 



Pa(a,,)-y-e-'^-- , (29) 

for which the problem can be solved exactly. Together 
with the constraint fj = F, we obtain a matrix of 
m rows and n columns. By demanding that all fj > 
0, one can in principle follow the same analysis as for 
real packings; we then average over all possible random 
matrices and consider only solutions with all fj > 0. In 
appendix IbI we derive that, in the limit that n,m —> oo 
with a fixed ratio p — m/n, the distribution becomes 

P(/) = dp) e-^^-P^f e-b(p)/= ^ (30) 

where b{p) is an almost linear function that has 6(0) = 
and 5(1) = I/tt. For square matrices, i.e. p = 1, we thus 
find that P{f) is a pure Gaussian centered around / = 0. 
This is illustrated in Fig. 113b : to calculate the P{f) for 
these non-square matrices we have evaluated Eqs. (jB3p 
and (|B5I) by Monte Carlo simulation. Tuning p to zero, 
the pressure constraint is dominant and we retrieve the 
pure exponential behavior that is also discussed in ap- 
pendix ^ 

So, we find that the tail of P{f) is a mixture of a 
Gaussian and an exponential, depending on the aspect 
ratio p of the matrix. However, for any value of p it is 
monotonically decreasing, and we never observe the peak 
that is extremely robust for real packing matrices. 

A relevant question of course is whether a Gaussian 
distribution of all matrix elements is representative for a 
matrix that is based on a real system of particles. Such 
a "real" packing matrix is not only sparse but also has 
ttij G [^1,1] in such a way that Newton's third law is re- 
spected. Unfortunately, it becomes very hard to work out 
the integrations if Pa (fly) is not Gaussian (4Qj or when 
correlations between matrix elements are imposed. For 
those systems, we have to rely on numerical simulations. 

2. Numerical simulations 

To numerically sample the ensemble discussed above, 
one first has to average over a representative number of 
allowed / for each matrix A, and then repeat this for 
many different matrices. However, only very few of the 
generated matrices actually have solutions for which all 
fj > 0. We have therefore focused our study on square 
random matrices, for which the phase space consists of a 
single point and the numerical effort is thus reduced to 
inverting each matrix. Starting from a random matrix for 
which all fj > 0, we apply a Monte Carlo simulation pro- 
cedure in which attempts are made to change a randomly 
selected element of A (except the elements corresponding 
to the pressure constraint). Such attempts are accepted 
with a probability given by the conventional Metropolis 
acceptance/rejection rules |4l| . In this way, we are able 
to explore the phase space of random matrices for which 
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FIG. 13: (a) Numerical evaluation of P{f) for matrices with 
dimensions ranging from 100 x 2 (p 0) to 100 x 100 (p = 1) 
illustrating crossover from exponential to Gaussian behavior 
(compare to Eg. I3UII . (b) Force distributions obtained with 
n X n Gaussian random matrices (with pressure constraint) 
for different values of n (curves) . For n = 50 the force distri- 
bution obtained with matrix elements from a uniform distri- 
bution is included for comparison (crosses). 



FIG. 14: (a) Force distributions obtained with 30 x 30 random 
matrices (with pressure constraint), with increasing sparse- 
ness. The distributions for Iz = 30 and Iz — 20 are indis- 
tinguishable, but for smaller Iz we see that the distribution 
becomes broader, (b) Here we show, for fixed Iz = 5, the 
emergence of a power-law in P{f) for large, sparse n x n ma- 
trices. 



all fj ^ 0, for any distribution of the matrix elements 

Gij. 

It is important to note that this numerical procedure 
is not precisely equivalent to the analysis of the Gaussian 
random matrices presented above. The reason for this is 
that a flat or uniform measure is not uniquely defined for 
continuous variables: a nonlinear transformation of vari- 
ables gives rise to a Jacobian that affects this flat phase 
space density. Since the coupling between and fj is in- 
deed nonlinear, the flat measure is ambiguous. However, 
one can show that the measure of the numerical scheme 
differs by a factor det(^) from P (^f, Aj of Eq. HB1|) . and 

we have verified that including this "weight factor" in the 
simulations only mildly alters the P{f) for small matrices 
{n < 5) and practically disappears for larger matrices. 

Square random matrices ~ Let us start the discussion 
with n X n square random matrices like the ones used 
for the analytical calculation above. This means one of 
the rows of the matrix represents the pressure constraint 
and the others are taken from a Gaussian distribution. In 
the limit n — > oo these were shown to give rise to a (half) 
Gaussian force distribution, see Eq. with p = I- The 
numerical results for n = 5, 10, 50 are shown in Fig. 113b. 
The distribution for n = 50 is indeed a Gaussian, as ex- 
pected for n ^ oo. The case n = 5 displays a very small 
peak at finite /, but this effect disappears quickly when 
n increases. Furthermore, Fig. 113b shows that the dis- 
tribution obtained with Gaussian matrix elements only 
slightly differs from the case of matrix elements taken 
from a flat distribution between — I and 1. 

Sparse matrices - A property of real packing matrices 
that is not represented by the random matrices is their 
sparseness: only those forces that push directly onto a 
given particle contribute to the force balance, and hence, 
most matrix elements are zero. On average, each row con- 
tains z nonzero elements, where z denotes the average co- 
ordination number. In order to investigate whether this 
sparseness is responsible for the non-monotonic P{f), we 
have generated a simple class of sparse random matri- 



ces: The matrices used are again n x n, but now with 
only Iz nonzero (Gaussian) elements per row (again, we 
leave the elements of the pressure contraint unaltered). 
These nonzero elements are arranged in a band-matrix 
like form. 

Force distributions for n = 30 and several values of 
Iz can be seen in Fig. I14h .. The maximum value of the 
distributions remains at / = and, surprisingly, it even 
increases as the matrix is more sparse. Uniformly dis- 
tributed elements gave almost identical results. It thus 
appears that the characteristic peak of P{f) is not di- 
rectly related to the sparseness of the matrix. In addition 
we found that for large sparse matrices, the tail of P{f) 
develops power-law scaling (Fig.lTlb). 

This demonstrates that a wide range of force distribu- 
tions can be obtained by varying the matrix properties, 
and that there is no simple answer to the question what 
properties of the matrix A are necessary to mimic realis- 
tic packings. In the light of this discussion, let us make 
the following remark. Recently, Ngan obtained a va- 
riety of force distributions similar to those obtained for 
real packing matrices in Sec. IIIBI and compatible with 
the form of Eq. (|30|l . These have been derived by mini- 
mizing an entropy functional under a pressure constraint 
similar to the one used in this paper |42j |. but without 
specifying the local microscopic equations of force bal- 
ance. One may therefore wonder whether it is possible 
to make a connection between the force ensemble and 
Ngan's work. On the other hand, the results of this sec- 
tion clearly illustrate that properties of the local equa- 
tions, which are absent in Ref. [25j, do play a crucial role: 
it can change P{f) from Gaussian to power-law. 

B. Perturbing a physical packing matrix 

In the previous section, we have shown that introduc- 
ing elements from real packing matrices to random square 
matrices does not easily lead to the characteristic peak 
in the resulting P{f). Therefore, we now investigate the 
reverse route, i.e. perturbing a real packing matrix by 
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FIG. 15: Definition of a rattler. Tfie net force on this rattler 
can only be zero if all forces involving this particle are zero. 
This means that the maximum angle between bonds, f3, is 
larger than tt. Such rattlers can arise when bonds are deleted 
or when the contact angles are randomly rotated (see text). 



slowly introducing randomness in the matrix elements. 
We perform three sorts of perturbations. In the first, 
the angles of the contacts are randomly varied, which 
ensures that the topology of the contact network remains 
unaltered. In the second, we randomly delete contacts, in 
the third, we randomly add contacts. In all three cases, 
the P{f) loses its maximum for sufficiently strong per- 
turbation. We show how for the first two protocols, this 
behavior appears to be correlated to the emergence of 
"rattlers" (see Fig. [T3|) . 

We have first constructed a matrix corresponding to an 
irregular packing of 1024 bi-disperse disks (50:50 mixture, 
size ratio 1.4) by molecular dynamics simulations using 
a 12-6 Lennard- Jones potential with the attractive tail 
cut off Til • This system is quenched below the glass 
transition {ksTg ~ 1.1) and its energy is minimized using 
a steepest descent algorithm, which guarantees that there 
is at least one stable force network. The resulting packing 
consists of 2814 bonds so z ~ 5.5. 

The effects on the P{f) for the force ensembles cor- 
responding to the perturbed matrices is illustrated in 
Fig. ^1 In Fig. I15b-b we illustrate the effect of rotat- 
ing the contacts by random angles uniformly generated 
between — and Ac/). With increasing A(/), the packing 
is getting more and more unphysical (corresponding less 
and less to a system of non-overlapping particles). Nev- 
ertheless, the topology of the network always remains 
the same, and Newton's third law is always respected. 
The resulting force distributions are computed using the 
algorithm described in Sec. HTIand averaged over all ran- 
domly generated perturbations of our original matrix A. 
In Fig. 116b we have plotted P{f) for different values of 
A(/). For small A0 we obtain the characteristic shape 
of P{f) for jammed systems similar to Fig. Small 
perturbations (A(/) < 0.2 rad) hardly change P{f), but 
at larger A(/) the peak around (/) disappears and P{f) 
looks "unjammed". For A0 > 0.75 we were no longer 
able to obtain solutions with all fj>0. 

This clearly shows that the conditions of a sparse 
matrix respecting the packing topology, elements dis- 
tributed between [—1, 1], and the incorporation of New- 
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FIG. 16: Variation of P{f) and number of rattlers when per- 
turbing a realistic packing matrix, (a-b) Variations of the con- 
tact angle randomly selected from [— A0, A0]; P{f) evolves 
from peaked to monotonic (a). The density of rattlers p (open 
symbols), and the rms variation of P{f) (stars) with respect 
to the unperturbed situation are roughly proportional (b). A 
similar scenario occurs when bonds are randomly deleted (c- 
d). When bonds are added, however, no rattlers are created 
but P{f) still evolves to a monotonic form (e-f). 



ton's third law into A are not sufficient to obtain the 
characteristic peak in P{f )- Even at relatively small per- 
turbations of A the shape of P{f) changes quite abruptly. 
Furthermore, our simulations clearly show that we are 
not even guaranteed to find a solution of the problem for 
a randomized matrix: only a very small fraction of all 
possible matrices lead to a solution for which all fj > 0. 
So, even though the emergence of a non-monotonic P{f) 
is extremely robust for packing matrices, it appears to 
be not at all a generic feature for arbitrary matrices. 

The amount of rattlers (Fig. ll5|) due to the randomiza- 
tion of the angles is small, but can be seen as a crude mea- 
sure of the contact geometry. To our suprise, the evolu- 
tion of the average amount of rattlers, and the rms devia- 
tion of P{f) from the unperturbed distribution arc fairly 
proportional (Fig. 116b ). Here, this rms deviation has 



been measured as y / df{Po{f) — P(/))2, where Po{f) 
denotes the unperturbed distribution. 

When bonds are deleted, a similar scenario occurs. 
Again the P(/)'s loose their peak and the rms devia- 
tion of P{f) follows the amount of rattlers quite well 
(Fig. ITBb-d) . On the other hand, when bonds are added, 
no rattlers are generated, but the P{f) still exhibit the 
same trend (Fig.llGb-f). Curiously, all the P(/)'s for the 
cases of added contacts appear to intersect in two points 



14 



lO'^ 
10-1 



S 



10" 

10" 

10" 

10" 
10-' 





(a) _ 




6 balls 






- - 15 balls 






- - 55 balls 






210 balls 







10« 
10-1 
10-2 
10-3 
10-4 

10-^ 
10-6 





(b) 


6 >v 




- - 15 \ 




--- 55 




— 210 





1 2 



4 5 6 



5 



10 15 20 25 30 35 



FIG. 17: Logarithmic plots of the P{f) for snooker triangles 
of increasing size as function of / (a) and (b) illustrate that 
the tails of these distribution decay faster than exponential 
but slower than Gaussian. 



fFig. 116( 3): we have no explanation for this phenomenon. 



VI. DISCUSSION 

In this paper we have proposed a novel ensemble ap- 
proach to athermal hard particle systems. The full set 
of mechanical equilibrium constraints were incorporated, 
in contrast to more local approximations or force chain 
models jM, 15, J5.,_26, 43, 44J. The basic idea is to exploit 
the separation of force and packing scales by simply aver- 
aging with equal probability over all mechanically stable 
force configurations for a fixed contact geometry. There 
are thus two important ingredients, namely the assump- 
tion of a flat (Edwards-like) measure in the force space 
and the fact that packings are hyperstatic. As the flat 
or uniform measure can not be justified from first prin- 
ciples, the emerging force probability distribution P{f) 
provides a first important test. For small forces, the 
ensemble nicely reproduces the typical non-monotonic 
behavior that has been found in numerous experiments 
and numerical simulations. Also, P{f) remains finite at 
f = 0, which has been the problem of earlier models 

lHlii. 

Let us now discuss the tails of the distribution. From 
Eq. Q we can only predict the asymptotic behavior of 
the slowest decaying term, corresponding to the minimal 
value of h\. For 2D packings one can show that this min- 
imal b\ cx 1/ViV, but since D (x N, the contribution to 

P{f) of this term decays as e~^^; this term thus pro- 
vides a sharp cutoff close to the maximal force. However, 
there will be a distribution of 6a 's, and in order to resolve 
the tail of P{f) one would really have to know all coeffi- 
cients in Eq. ((Jj) for large enough systems. In Fig. 1171 we 
again plot the numerically obtained P{f) for snooker- 
packings. Although the systems are of limited size, it 
appears that the distributions have tails that neither are 
purely exponential nor purely Gaussian. The differences 
are subtle, and may be sensitive to numerical details. 
Even though numerical and analytical distributions for 
small packings appear to be in perfect agreement on a lin- 
ear scale, on similar log scales the numerical curve seems 
to slightly underestimate the large fluctuations. While 



the numerical precision is about 10~^ around (/), the 
relative differences between numerical and exact results 
become about 5% around / = 4(/). In the literature, 
there has recently been some debate on the true nature 
of the tails 45] : while the carbon paper experiments un- 
doubtedly yield exponential tails, it appears that most 
numerically obtained P{f) display some downward cur- 
vature when plotted on log-lin scales. It has also been ar- 
gued that individual packings are not self-averaging and 
that tails appear Gaussian or e xpon ential depending on 
how the ensemble is normalized At present, we can 
therefore neither confirm nor falsify the validity of the 
fiat measure based on the tail of 

Unger et al. 0| recently proposed another test for 
the uniform measure. For frictional packings, they com- 
pared force configurations that emerge in a dynamical 
process to those obtained from a random sampling of the 
force space. They found that the dynamical solutions 
are located more centrally within the force space, and 
therefore concluded that the fiat measure does not apply. 
While this is definitely an interesting observation, this 
claim strongly depends on the "flatness" of their numeri- 
cal sampling of the solution space, for which no evidence 
is provided. Counterintuitively, if the physical force net- 
works were indeed more central, the ensemble P{f) would 
even overestimate the large force fluctuations. Therefore, 
the validity of the flat measure remains an open issue. 

A second important ingredient of the force network 
ensemble is that there is no unique force solution for 
a given contact network, i.e. packings are hyperstatic. 
While most packings are indeed hyperstatic, the precise 
degree of indeterminacy may depend on material param- 
eters and construction history |l9t 0| . It appears 
that strict isostaticity is only found for infinitely hard 
particles without friction, or with unphysically large fric- 
tion coefficients. The present study was performed with 
highly coordinated regular packings, which are more hy- 
perstatic than most physical packings. The coordination 
number is therefore an important parameter that remains 
to be explored. It may very well be that the predictive 
power of the ensemble depends on this degree of indeter- 
minacy. 

Metabolic networks - While preparing this paper, we 
have become aware of a striking analogy between the 
force ensemble and the problem of metabolic networks 
I27L |2^ . These are networks of biochemical reactions, 
in which the metabolite concentrations (particle posi- 
tions) and the reaction and transport fluxes (interpar- 
ticle forces) are the variables of the problem. In princi- 
ple the fluxes follow from the concentrations, similar to 
how the forces follow from the particle positions. This 
coupling involves intricate reaction-diffusion dynamics, 
for which numerical values of most rates are not known. 
In practice, however, a separation of time-scales occurs: 
the metabolite concentrations quickly adjust (seconds) 
to global changes in the network (minutes) jig. Very 
much like we employed the separation of length-scales, a 
successful strategy has been to treat the fluxes as inde- 
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pendent variables and resolve the steady state from the 
stoichiometry of the network. 

Mathematically, the problem then reduces to an under- 
determined matrix problem with non-negative flux vari- 
ables, which is identical to the equations defining the 
force network ensemble. It turns out that for different 
metabolic maps the number of fluxes is always larger 
than the number of metabolites and therefore these sys- 
tems are "hyperstatic" . The main difference with respect 
to the force problem, however, is the network structure 
defining the matrix: metabolic networks are scale-free, 
i.e. with highly uneven connectivities. This leads to 
power-law flux distributions , which is very different 
from the P{f) within the force ensemble. This touches 
upon the deep question of how network statistics relate 
to the underlying network structures. In Sec. we have 
found that, indeed, P{f) can range from Gaussian to 
power-law when changing the properties of the matrix 
defining the ensemble. 

For metabolic networks the main interest is to find 
solutions in which the production of "biomass" is opti- 
mized. In contrast to the averaging procedure within 
the force ensemble, this corresponds to finding the "ex- 
treme pathways" that form the corners of the hyperpoly- 
gon defining the solution space 27]. In fact, the force 
network solutions shown in Figs. 13 and 1911 II are such 
extreme pathways. We speculate that a systematic anal- 
ysis of extreme solutions may give additional insight in 
the geometrical properties of the phase space and the 
emergent force statistics - see also Ref. |2J|. It would 
furthermore be interesting to see whether for disordered 
packings there still exist localized linearly propagating 
solutions such as shown in Fig. O or whether all particles 
have to be involved into the force network. 

Outlook - A number of crucial questions can possi- 
bly be addressed within our framework. (1) By separat- 
ing the contact geometry from the forces, we can start 
to disentangle the separate roles of contact and stress 
anisotropics in sheared systems. In particular, we have 
already shown that the ensemble comprises an "unjam- 
ming" transition for shear stresses above a critical value 
[Tl|. Furthermore, the contact and force networks ex- 
hibit different anisotropics under different construction 
histories 0, 0|. We suggest that the contact network 
anisotropics may be sufficient to obtain the pressure dip 
under sand piles. (2) As also illustrated by |23,|2J|, our 
approach is perfectly suited to include frictional forces. 
While these forces are difficult to express in a force law, 
they are easy to constrain by the Coulomb inequality. 
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APPENDIX A: THE PRESSURE CONSTRAINT 

In this appendix we first show that the sum of all forces 
fij is constant for regular packings under a fixed exter- 
nal pressure. This provides a conservation law similar to 
the conservation of total energy in the microcanonical en- 
semble. We therefore revisit the problem of partitioning 
a conserved quantity in the second part of this appendix. 
One can compute the stress from the contact forces fij 

as 



(Al) 



where V represents the volume of the system 47]. The 
vector Tij = — rj denotes the interparticle distance, 
which for monodisperse particles of diameter d always 
has |r| = d. For packings of frictionless disks, one can 
therefore write 



— \ ^ f 2 

^xx — / ^ Jij cos (pij , 



{y} 



(A2) 
(A3) 



{y} 



where ifij indicates the angle of the contact with repect 
to the horizonal x-axis. Taking the trace of the stress 

tensor, we find axx + <^yy ~ ^Y fij- So indeed, a con- 
stant pressure condition is equivalent to a constraint for 
the sum of all contact forces, at least for monodisperse 
packings. To a good approximation this remains valid 
for polydisperse packings, since in practice, the forces 
are uncorrelated to |rl so that (|r|) can be taken out of 
the sum in Eq. [tIITH. 

Let us now consider the statistical properties of a set 
of n independent non-negative variables Xj > 0, that is 
contrained by a conservation law 



(A4) 



The phase space of these variables is a (n—l) dimensional 
simplex of volume 







n 

(n- 1)! ' 



dxi 




(A5) 



where the integrals can be evaluated e.g. by Fourier rep- 
resentation of the (5-function. Assigning an equal proba- 
bility to all sets {xj} obeying Eq. ljA4p . we compute the 
probability of a single variable P{x; n) as 
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P{x; n) 



r!„(X) 



J=2 



Q{X ~ x) 



(A6) 



In the thermodynamic limit n — > oo, this becomes the 
purely exponential "Boltzmann" distribution 



P(x;oo) = T^e' 
a; 



-x/{x) 



(A7) 



where {x) — X/n. For finite n > 2, however, this distri- 
bution is always monotonically decreasing. 

In this paper we encounter two (small) packing con- 
figurations for which the force network ensemble can be 
reduced to the simple problem discussed above, so that 
force distributions of the type Eq. (|A6ll are found - see 
Figs. 0] and ^| In general, however, the constraints of 
force balance on each particle are more complicated and 
lead to non-monotonic P{f). 



APPENDIX B: DERIVATION OF THE 
GAUSSIAN RANDOM MATRIX P(/) 

In this appendix we show how Eq. (|3U|1 is obtained. 
We study the ensemble defined by 

m I n \ 

xIlM^^^^'^M ' ^^^^ 

i=2 \j=l j 

where we define 

In order to be consistent with the notation in Sec. lIIII we 
reserve the index i ~ I for the inhomogeneous pressure 
constraint. The force distribution P{fj) becomes 



(B3) 



where 



da.; 



(B4) 



The advantage of taking Gaussian elements aij is that 
they can be integrated out explicitly, using the Fourier 
representations of Eq. H13|l : 



/m / ^ 
dAPa{A)X{5[Y^a,,f, 

1=2 = l 

rOQ 1 poo 

n/ frll/ 



iSi O-ij fj 



(m-l)/2 



J ■'3 



(B5) 



It is convenient to bring the factor /J to the exponent 
using the relation 



1 1 



dtr~^e 



T{k) 7o 

Introducing this auxilary variable t, P{fj) becomes 



(B6) 



X / dtt^ e 











where we have used a shorthand 

POO 



(B7) 



(B8) 



We now exponentiate the full integrand of Eq. ljB7|l . so 
that P{fj) becomes 

/OO 7 />OC 
^ / dte--^^-*/^*(--*) , (B9) 
-OO Jo 



where 



$(iisi,t) = isiF+ ( ^^^^—^ ] ln(t) 



+ l)ln[g(zsi,t)] . (BIO) 

If we now fix (/) = 1 by taking F — n, one observes 
that all terms of the phase $ are extensive in n or m. 
In the limit where both n, m — > oo, we can thus evalu- 
ate the integrals using a saddle-point approximation. By 
determining the stationary phase, i.e. d^/dsi = and 
d^/dt — 0, one finally arrives at the result of Eq. H30I) : 



P(/) = c(p)e-(i-'')^e-''('')^' 



(Bll) 



The function b{p) varies almost linearly between 6(0) = 
and 6(1) = I/tt. 
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